Hive Parameter Analysis

library(MASS)
library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6     ✔ purrr   0.3.4
✔ tibble  3.1.7     ✔ dplyr   1.0.9
✔ tidyr   1.2.0     ✔ stringr 1.4.0
✔ readr   2.1.2     ✔ forcats 0.5.1
── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
✖ dplyr::select() masks MASS::select()

library(ggplot2)
library(readr)
library(viridis)

Loading required package: viridisLite

library(ggfortify)
library(bbmle) #For ICtab

Loading required package: stats4

Attaching package: ‘bbmle’

The following object is masked from ‘package:dplyr’:

    slice

library(car)

Loading required package: carData

Attaching package: ‘car’

The following object is masked from ‘package:dplyr’:

    recode

The following object is masked from ‘package:purrr’:

    some

library(emmeans)

 
round_any <- function(x, accuracy, f=round){f(x/ accuracy) * accuracy}

round_any <- function(x, accuracy, f=round){f(x/ accuracy) * accuracy}

# trial_types <- c("none","brood","worker","queen","all")
# 
# for (t in trial_types){
#   bee_param_df <- list.files(path=paste0("hive_data/heat_",t,"/"), full.names = TRUE) %>%
#                 lapply(read_csv, show_col_types = FALSE) %>%
#                 bind_rows %>%
#                 mutate(filenames = list.files(path=paste0("hive_data/heat_",t,"/"), full.names = TRUE))
# 
#   write.csv(bee_param_df,paste0("bee_data_isaac_heat_",t,".csv"),row.names = FALSE)
# }

What these all mean: n = queen cells per hour rb = brood radius rn = necter radius w = total daily honey pph = pollen ratio ph = honey consumption ratio pp = pollen consumption ratio k = consumption probability value

broodMetric = average number of brood surrounding other brood pollenRing = average min distance between honey and brood

brood_bee_param <- read.csv("bee_data_isaac_heat_brood.csv") %>% mutate(type = "Brood") %>% mutate(wha = 0, qha = 0)
worker_bee_param <- read.csv("bee_data_isaac_heat_worker.csv") %>% mutate(type = "Worker") %>% mutate(bhd = 0, qha = 0)
queen_bee_param <- read.csv("bee_data_isaac_heat_queen.csv") %>% mutate(type = "Queen") %>% mutate(bhd = 0, wha = 0)
none_bee_param <- read.csv("bee_data_isaac_heat_none.csv") %>% mutate(type = "None") %>% mutate(bhd = 0,wha = 0, qha = 0)
all_bee_param <- read.csv("bee_data_isaac_heat_all.csv") %>% mutate(type = "All")

bee_heat_param_df <- rbind(brood_bee_param,worker_bee_param) %>% rbind(.,queen_bee_param) %>% rbind(.,all_bee_param) %>% rbind(.,none_bee_param)

ggplot(bee_heat_param_df, aes(x = type, y = pBroodHeat))+
  geom_boxplot()+
  theme_classic()

heat_glm <- glm(pBroodHeat ~ 0 + type, data = bee_heat_param_df)

contrast(emmeans(heat_glm, "type"), "pairwise", adjust = "Tukey")

 contrast       estimate      SE   df t.ratio p.value
 All - Brood     -0.0545 0.00819 3835  -6.650  <.0001
 All - None      -0.4391 0.00819 3835 -53.600  <.0001
 All - Queen     -0.2398 0.00819 3835 -29.271  <.0001
 All - Worker    -0.3607 0.00819 3835 -44.021  <.0001
 Brood - None    -0.3847 0.00772 3835 -49.798  <.0001
 Brood - Queen   -0.1853 0.00772 3835 -23.993  <.0001
 Brood - Worker  -0.3062 0.00772 3835 -39.638  <.0001
 None - Queen     0.1993 0.00772 3835  25.805  <.0001
 None - Worker    0.0785 0.00772 3835  10.160  <.0001
 Queen - Worker  -0.1208 0.00772 3835 -15.645  <.0001

P value adjustment: tukey method for comparing a family of 5 estimates 

ggplot(all_bee_param %>% na.omit(), aes(x = qha, y = bhd, color = pBroodHeat))+
  geom_point()+
  scale_color_viridis()+
  theme_classic()

rounded_all_bee_param <- all_bee_param %>% mutate(rounded_qha = round_any(qha, 0.25),
                                                 rounded_bhd = round_any(bhd, 0.25))

ggplot(rounded_all_bee_param %>% na.omit(), aes(x = rounded_qha, y = rounded_bhd, z = pBroodHeat))+
  geom_contour_filled()+
  theme_classic()

NA
NA

Looking at correlation between heat parameter and brood in the place

ggplot(bee_heat_param_df %>% filter(bhd > 0), aes(x = bhd, y = pBroodHeat))+
  geom_point()+
  geom_smooth(method = 'glm', formula = 'y ~ x')+
  theme_classic()

cor.test((bee_heat_param_df %>% filter(bhd > 0))$bhd,(bee_heat_param_df %>% filter(bhd > 0))$pBroodHeat)

    Pearson's product-moment correlation

data:  (bee_heat_param_df %>% filter(bhd > 0))$bhd and (bee_heat_param_df %>% filter(bhd > 0))$pBroodHeat
t = -32.341, df = 1438, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.6778413 -0.6179628
sample estimates:
       cor 
-0.6489056 

ggplot(bee_heat_param_df %>% filter(qha > 0), aes(x = qha, y = pBroodHeat))+
  geom_point()+
  geom_smooth(method = 'glm', formula = 'y ~ x')+
  theme_classic()

cor.test((bee_heat_param_df %>% filter(qha > 0))$qha,(bee_heat_param_df %>% filter(qha > 0))$pBroodHeat)

    Pearson's product-moment correlation

data:  (bee_heat_param_df %>% filter(qha > 0))$qha and (bee_heat_param_df %>% filter(qha > 0))$pBroodHeat
t = -27.375, df = 1438, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.6182844 -0.5503031
sample estimates:
       cor 
-0.5853215 

ggplot(bee_heat_param_df %>% filter(wha > 0), aes(x = wha, y = pBroodHeat))+
  geom_point()+
  geom_smooth(method = 'glm', formula = 'y ~ x')+
  theme_classic()

cor.test((bee_heat_param_df %>% filter(wha > 0))$wha,(bee_heat_param_df %>% filter(wha > 0))$pBroodHeat)

    Pearson's product-moment correlation

data:  (bee_heat_param_df %>% filter(wha > 0))$wha and (bee_heat_param_df %>% filter(wha > 0))$pBroodHeat
t = -5.9184, df = 1438, p-value = 4.058e-09
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2042354 -0.1033708
sample estimates:
       cor 
-0.1542049 

Let’s make some stepwise models! Just with the all data for now so there aren’t so many zeroes

hive_pca <- prcomp(bee_heat_param_df %>% select(-c(...1, trial_n, type, days, n,
                                                   pBrood, pPollen, pHoney, pEmpty,
                                                   pBroodHeat, pPollenHeat, pHoneyHeat, pEmptyHeat,
                                                   broodMetric,pollenRing)),
                   center = TRUE,scale. = TRUE)

Error in colMeans(x, na.rm = TRUE) : 'x' must be numeric

heat_brood_glm <- glm(pBroodHeat ~ ., data = all_bee_param %>% dplyr::select(-c(...1, trial_n, type, days, n,
                                                   pBrood, pPollen, pHoney, pEmpty,
                                                   pPollenHeat, pHoneyHeat, pEmptyHeat,
                                                   broodMetric,pollenRing)))

stepped_model <- stepAIC(heat_brood_glm, direction = "both")

Start:  AIC=-1644.08
pBroodHeat ~ rb + rn + w + pph + ph + pp + k + bhd + wha + qha

       Df Deviance     AIC
- wha   1   2.7657 -1646.0
- rb    1   2.7660 -1646.0
- k     1   2.7710 -1644.8
- pph   1   2.7733 -1644.3
<none>      2.7655 -1644.1
- rn    1   2.7757 -1643.7
- w     1   2.7808 -1642.5
- ph    1   2.7828 -1642.1
- pp    1   2.8025 -1637.6
- qha   1   3.6132 -1475.0
- bhd   1   4.1686 -1383.5

Step:  AIC=-1646.02
pBroodHeat ~ rb + rn + w + pph + ph + pp + k + bhd + qha

       Df Deviance     AIC
- rb    1   2.7662 -1647.9
- k     1   2.7711 -1646.8
- pph   1   2.7739 -1646.2
<none>      2.7657 -1646.0
- rn    1   2.7759 -1645.7
- w     1   2.7812 -1644.5
+ wha   1   2.7655 -1644.1
- ph    1   2.7832 -1644.0
- pp    1   2.8028 -1639.5
- qha   1   3.6153 -1476.6
- bhd   1   4.1692 -1385.3

Step:  AIC=-1647.9
pBroodHeat ~ rn + w + pph + ph + pp + k + bhd + qha

       Df Deviance     AIC
- k     1   2.7714 -1648.7
- pph   1   2.7743 -1648.0
<none>      2.7662 -1647.9
- rn    1   2.7764 -1647.5
- w     1   2.7821 -1646.2
+ rb    1   2.7657 -1646.0
- ph    1   2.7833 -1646.0
+ wha   1   2.7660 -1646.0
- pp    1   2.8029 -1641.5
- qha   1   3.6153 -1478.6
- bhd   1   4.1692 -1387.3

Step:  AIC=-1648.71
pBroodHeat ~ rn + w + pph + ph + pp + bhd + qha

       Df Deviance     AIC
<none>      2.7714 -1648.7
- pph   1   2.7802 -1648.7
- rn    1   2.7807 -1648.6
+ k     1   2.7662 -1647.9
- ph    1   2.7879 -1646.9
+ rb    1   2.7711 -1646.8
+ wha   1   2.7713 -1646.7
- w     1   2.7888 -1646.7
- pp    1   2.8093 -1642.0
- qha   1   3.6258 -1478.7
- bhd   1   4.1708 -1389.1

Anova(stepped_model)

Analysis of Deviance Table (Type II tests)

Response: pBroodHeat
    LR Chisq Df Pr(>Chisq)    
rn      2.11  1   0.146522    
w       3.96  1   0.046605 *  
pph     2.00  1   0.157613    
ph      3.77  1   0.052131 .  
pp      8.63  1   0.003303 ** 
bhd   319.12  1  < 2.2e-16 ***
qha   194.83  1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ggplot(all_bee_param, aes(bhd, pBroodHeat, color = pBroodHeat))+
  geom_point()+
  scale_color_viridis()+
  theme_classic()

ggplot(all_bee_param, aes(qha, pBroodHeat, color = pBroodHeat))+
  geom_point()+
  scale_color_viridis()+
  theme_classic()

now let’s see about the order of those models

bhd_m <- glm(pBroodHeat ~ bhd, data = all_bee_param)
bhd_2_m <- glm(pBroodHeat ~ poly(bhd,2), data = all_bee_param)
bhd_3_m <- glm(pBroodHeat ~ poly(bhd,3), data = all_bee_param)
bhd_ln_m <- glm(pBroodHeat ~ log(bhd), data = all_bee_param)

ICtab(bhd_m,bhd_2_m,bhd_3_m,bhd_ln_m)

         dAIC df
bhd_2_m   0.0 4 
bhd_3_m   1.3 5 
bhd_ln_m 15.3 3 
bhd_m    43.6 3 

ggplot(all_bee_param, aes(bhd, pBroodHeat, color = pBroodHeat))+
  geom_point()+
  geom_smooth(method = "glm", formula = "y ~ poly(x,2)", se = FALSE) +
  scale_color_viridis()+
  theme_classic()

qha_m <- glm(pBroodHeat ~ qha, data = all_bee_param)
qha_2_m <- glm(pBroodHeat ~ poly(qha,2), data = all_bee_param)
qha_3_m <- glm(pBroodHeat ~ poly(qha,3), data = all_bee_param)
qha_ln_m <- glm(pBroodHeat ~ log(qha), data = all_bee_param)

ICtab(qha_m,qha_2_m,qha_3_m,qha_ln_m)

         dAIC df
qha_m     0.0 3 
qha_3_m   0.3 5 
qha_2_m   1.9 4 
qha_ln_m 31.9 3 

ggplot(all_bee_param, aes(qha, pBroodHeat, color = pBroodHeat))+
  geom_point()+
  geom_smooth(method = "glm", formula = "y ~ x", se = FALSE) +
  scale_color_viridis()+
  theme_classic()

From the none data, looking at the good ones. The pollen ring metric seems a bit off

Step AIC to see what predicts good or bad

Anova(stepped_good_model)

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Analysis of Deviance Table (Type II tests)

Response: good_or_bad
   LR Chisq Df Pr(>Chisq)    
rn  22.3956  1  2.219e-06 ***
k    7.7241  1   0.005449 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ggplot(none_bee_param_good_bad, aes(x = n, fill = good_or_bad))+
  geom_histogram(bins = 10)+
  theme_classic()

ggplot(none_bee_param_good_bad, aes(x = rb, fill = good_or_bad))+
  geom_histogram(bins = 10)+
  theme_classic()

ggplot(none_bee_param_good_bad, aes(x = rn, fill = good_or_bad))+
  geom_histogram(bins = 10)+
  theme_classic()

ggplot(none_bee_param_good_bad, aes(x = w, fill = good_or_bad))+
  geom_histogram(bins = 10)+
  theme_classic()

ggplot(none_bee_param_good_bad, aes(x = pph, fill = good_or_bad))+
  geom_histogram(bins = 10)+
  theme_classic()

ggplot(none_bee_param_good_bad, aes(x = ph, fill = good_or_bad))+
  geom_histogram(bins = 10)+
  theme_classic()

ggplot(none_bee_param_good_bad, aes(x = pp, fill = good_or_bad))+
  geom_histogram(bins = 10)+
  theme_classic()

ggplot(none_bee_param_good_bad, aes(x = k, fill = good_or_bad))+
  geom_histogram(bins = 10)+
  theme_classic()

OLDER CODE, looking at model type, we decided on M2 ####

rand_bee_param <- read.csv("bee_data_isaac_rand.csv") %>% mutate(type = "Rand")
m1_bee_param <- read.csv("bee_data_isaac_m1.csv") %>% mutate(type = "M1") %>% select(-trial_n)
m2_bee_param <- read.csv("bee_data_isaac_m2.csv") %>% mutate(type = "M2") %>% select(-trial_n)
m2E_bee_param <- read.csv("bee_data_isaac_m2_empty.csv") %>% mutate(type = "M2E") %>% select(-trial_n)
m2K_bee_param <- read.csv("bee_data_isaac_m2_katie.csv") %>% mutate(type = "M2K") %>% select(-trial_n)

bee_param_df <- rbind(m1_bee_param,m2_bee_param) %>% rbind(.,rand_bee_param) %>% rbind(.,m2E_bee_param) %>% rbind(.,m2K_bee_param)

bee_param_df_perc <- bee_param_df %>% select(c(type,pBrood,pHoney,pPollen,pEmpty)) %>%
                                      pivot_longer(c(pBrood,pHoney,pPollen,pEmpty))

ggplot(bee_param_df, aes(x = broodMetric, y = pollenRing, color = type))+
  geom_point(alpha = 0.25)+
  ylim(0,16)+
  xlim(0,6)+
  theme_classic()

ggplot(bee_param_df_perc, aes(x = type, y = value, fill = name))+
  geom_boxplot()+
  theme_classic()

ggplot(m1_bee_param, aes(x = pPollen, y = pph))+
  geom_point()+
  theme_classic()

broodMetricGLM_m0 <- glm(broodMetric ~ 1, data = m2_bee_param)
broodMetricGLM_m1 <- glm(broodMetric ~ n + rb + rn + w + pph + ph + pp + k, data = m2_bee_param)
broodMetricGLM_step <- step(broodMetricGLM_m0, direction = "both", scope = formula(broodMetricGLM_m1), trace = 0)

summary(broodMetricGLM_step)

pollenRingGLM_m0 <- glm(pollenRing ~ 1, data = m2_bee_param)
pollenRingGLM_m1 <- glm(pollenRing ~ n + rb + rn + w + pph + ph + pp + k, data = m2_bee_param)
pollenRingGLM_step <- step(pollenRingGLM_m0, direction = "both", scope = formula(pollenRingGLM_m1), trace = 0)

summary(pollenRingGLM_step)

queen_pos_df <- read.csv("queen_pos_data_28_5.csv")

queen_pos_df_heatmap <- queen_pos_df %>% mutate(QueenXRound = round_any(QueenX,5),
                                                QueenYRound = round_any(QueenY,5)) %>%
                                        group_by(QueenXRound,QueenYRound) %>%
                                        summarise(meanAngleMoved = mean(AngleMoved),
                                                  count = n())

`summarise()` has grouped output by 'QueenXRound'. You can override using the `.groups` argument.

ggplot(queen_pos_df, aes(x = QueenX, y = QueenY))+
  geom_path()+
  geom_point(aes(x=11,y=48.5),colour="red") +
  scale_color_viridis() +
  theme_classic()

ggplot(queen_pos_df, aes(x = QueenX, y = QueenY, color = Dist2Heat))+
  geom_point()+
  geom_point(aes(x=11,y=48.5),colour="red") +
  scale_color_viridis() +
  theme_classic()

ggplot(queen_pos_df, aes(x = QueenX, y = QueenY, color = cos(AngleOppHeat)))+
  geom_point()+
  geom_point(aes(x=11,y=48.5),colour="red") +
  scale_color_viridis() +
  theme_classic()

ggplot(queen_pos_df, aes(x = QueenX, y = QueenY, color = sin(AngleOppHeat)))+
  geom_point()+
  geom_point(aes(x=11,y=48.5),colour="red") +
  scale_color_viridis() +
  theme_classic()

ggplot(queen_pos_df_heatmap, aes(x = QueenXRound, y = QueenYRound, color = cos(meanAngleMoved)))+
  geom_point(size = 5)+
  geom_point(aes(x=11,y=48.5),colour="red") +
  scale_color_viridis() +
  theme_classic()

ggplot(queen_pos_df_heatmap, aes(x = QueenXRound, y = QueenYRound, color = sin(meanAngleMoved)))+
  geom_point(size = 5)+
  geom_point(aes(x=11,y=48.5),colour="red") +
  scale_color_viridis() +
  theme_classic()

ggplot(queen_pos_df_heatmap, aes(x = QueenXRound, y = QueenYRound, color = count))+
  geom_point(size = 5)+
  geom_point(aes(x=11,y=48.5),colour="red") +
  scale_color_viridis() +
  theme_classic()